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ON i 

ON • Abstract 

ON ■ 

t-H ' The coupled Einstein-Dirac equations for a static, spherically symmetric system 

of two fermions in a singlet spinor state are derived. Using numerical methods, we 
construct an infinite number of soliton-like solutions of these equations. The stability 
of the solutions is analyzed. For weak coupling (i.e., small rest mass of the fermions), 
all the solutions are linearly stable (with respect to spherically symmetric perturba- 
tions), whereas for stronger coupling, both stable and unstable solutions exist. For 
the physical interpretation, we discuss how the energy of the fermions and the (ADM) 
mass behave as functions of the rest mass of the fermions. Although gravitation 
Q> ' is not renormalizable, our solutions of the Einstein-Dirac equations are regular and 

t*^* . well-behaved even for strong coupling. 

^ ■ 1 Introduction 

ON ■ 

In recent years, there has been much interest in the coupling of Einstein's field equations 
to Yang- Mills equations. In this case, the attractive gravitational force is balanced by the 
^ | repulsive Yang-Mills force, and this interaction has led to many interesting and surprising 

results; see for example [|l]]-[|lO]]. In this paper, we consider the coupling of Einstein's 
t> ■ equations to the Dirac equation. Here the necessary repulsive mechanism is provided by 

^ , the Heisenberg Uncertainty Principle. 

$— i ' The Einstein-Dirac equations take the form 



C3 



Ri-lRS) = -87rT] , (G-m)tf = , (1.1) 

where Tj is the energy-momentum tensor of the Dirac particle, G denotes the Dirac op- 
erator (see |13[| ), and ^ is the wave function of a fermion of mass m. As in the above- 
mentioned earlier studies, we consider static, spherically symmetric solutions. Since the 
spin of a fermion has an intrinsic orientation in space, a system consisting of a single Dirac 
particle cannot be spherically symmetric. In order to maintain the spherical symmetry, 
we are led to the study of two fermions having opposite spin; i.e. to a singlet spinor state. 
Of course, such a configuration does not represent a realistic physical system due to the 
absence of the electromagnetic interaction. More precisely, neglecting the electromagnetic 
interaction corresponds to the limiting case where the masses of the fermions become so 
large (of the order (Planck length) -1 ) that the gravitational interaction becomes the dom- 
inant force. Nevertheless, we view this study as a model problem worth considering in 
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order to get some understanding of the equations and their solutions. In a future pub- 
lication, we will consider the more physically realistic situation where the Einstein-Dirac 
equations are coupled to an electromagnetic field (Maxwell's equations). 

Our results are based on a certain ansatz, whereby we reduce the 4-component Dirac 
spinors to a 2-component spinor system, <3? = (a, (3) with real functions a, j3. We show 
numerically that particle-like solutions of this type exist, both in the ground state, and 
in the excited states. For weak coupling, i.e. small mass m, the different solutions are 
characterized by the "rotation number", n = 0, 1,2, . . ., of the vector (a,f3) (we work in 
standard units h = c = Q = 1). The solution with n = is the ground state, and the 
solutions with n = 1,2,... describe the excited states. For small m, the solutions are 
(linearly) stable with respect to spherically symmetric perturbations. However, as m gets 
large, several states appear for each n. In fact, for every n, the mass spectrum (i.e., the 
plot of the binding energy vs. the rest mass) is a spiral curve which tends to a limiting 
configuration. This surprising result shows that for parameter values on this limiting 
configuration, there are an infinite number of excited states "in the n th mode" , while for 
parameter values near this limiting configuration, there are still a large, but finite number 
of such excited states. Furthermore, using topological methods and bifurcation theory (see 

part IV]), we show that in every mode, the stable solutions must become unstable as 
the binding energy increases. Although gravitation is not renormalizable (which means 
that the problem cannot be treated in a perturbation expansion), our solutions of the 
Einstein-Dirac equations are regular and well-behaved even for strong coupling. 



2 The Dirac Operator 

In this section, we shall derive the form of the Dirac operator in the presence of a static, 
spherically symmetric gravitational field. In preparation, we first give a brief mathemat- 
ical introduction of the Dirac theory in curved space-time. The Dirac operator G is a 
differential operator of first order 

G = i&W^j + B{x) , (2.1) 

where the Dirac matrices G J (x), (j = 0,1,2,3), and B(x) are (4 x 4) matrices, which 
depend on the space-time point x. The Dirac matrices and the Lorentzian metric are 
related by 

9 jk = \ G k ) , (2.2) 

where {., .} is the anti-commutator 

{G j , G k } = (G j G k + G k G j ) 

The basic difficulty with Dirac spinors in curved space-time is that, for a given Lorentzian 
metric, the Dirac matrices are not uniquely determined by the anti-commutation relations 
fl2.2| ). One way of fixing the Dirac matrices is provided by the formalism of spin and 
frame bundles (see e.g. the first section of flT| ). In this formulation, one chooses a frame 
(w a ) a= o,...,3 and represents the Dirac matrices as linear combinations of the Dirac matrices 
7 a of Minkowski space, 

&(x) = 5> J »7 a • 

a=0 
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The matrix B(x) is composed of the so-called spin connection coefficients, involving first 
partial derivatives of the metric and of the frame. It is quite common to choose for 
(u a ) a Newman- Penrose null frame; this choice is particularly convenient for metrics of 
Petrov type D (see for an introduction to the Newman-Penrose formalism and many 
applications, especially in the Kerr background). More generally, it is shown in [13| that 



all choices of Dirac matrices satisfying ( |2.2D yield unitarily equivalent Dirac operators. 
Furthermore, [13] gives explicit formulas for the matrix B in terms of the Dirac matrices 
G 3 . We prefer working with the formalism of [13] in the following, because it gives us 
more flexibility in choosing the Dirac matrices. 

The wave function ^ of a Dirac particle is a solution of the Dirac equation 

(G-m)* = . (2.3) 

On the wave functions, two different scalar products can be defined. In the first, we 
integrate the wave functions over all of space-time, 



<^ ]<£>> = J y$yj\g\drx , (2.4) 

where Vl/ = ^* ^ J ^ is the adjoint spinor (whose definition does not depend on the 

gravitational field; 0, 1 are (2 x 2) submatrices), and g denotes the determinant of the 
metric gjk- The scalar product ( |2.4[ ) is indefinite, but it is nevertheless useful to us because 
the Dirac operator is Hermitian with respect to it. The second scalar product is defined 
on the solutions of the Dirac equation. For this we choose a space-like hypersurface TL 
together with a (future-directed) normal vector field u, and set 

(tf|$)= / vi dp , (2.5) 

JH 

where dp is the invariant measure on the hypersurface H., induced by the metric gij. This 
scalar product is positive definite, and, as a consequence of the current conservation (cf. 
@) 

Vj^G^ = , (2.6) 

it is independent of the choice of the hypersurface 7i. In direct generalization of the 
expression tyj 0l fy in Minkowski space (see e.g |l|J), the integrand ^G 3 ^ Uj is interpreted 
as the probability density of the particle. This leads us to normalize solutions of the Dirac 
equation by requiring 

(tf | ) = 1 . (2.7) 



We now return to the Dirac operator ( |2.1| ). Suppose that a 4-dimensional space-time 
with metric gij is given. According to |l3|], we can choose for G 3 any 4x4 matrices which 
are Hermitian with respect to the scalar product Q2.4| ) and satisfy ( |2.2[ ). The matrix B(x) 
involves first derivatives of the Dirac matrices G 3 , and from [13], we have the explicit 
formulas 

B{x) = G 3 \x)E j {x) (2.8) 



with 



Ej = ~ p(d jP ) - Tr(G m V 3 G n ) G m G n + 1 Tr( P G, V m G m ) p (2.1 
Z lb o 



4! 



e ijkl G i G 3 G k G l (2.10) 
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(eijki is the totally antisymmetric tensor density). 

Now we will specify these formulas for the Dirac operator to static, spherically sym- 
metric space-times. In polar coordinates (t, r, ip), the metric can be written as (cf. 
©0) 

9ij = diag -i -r 2 , -r 2 sin 2 0) (2.11) 



1 1 



^ = diag (T 2 , -A ~, - ; ) (2.12) 

r z sin 



gf| = T _1 ,4~5 r 2 | sim?| 



with volume element 



where A = A{r) and T = T{r) are positive functions. We shall use this form of the metric 
to explicitly calculate the Dirac operator ( p.lj ). For the Dirac matrices G 3 (x), we take an 
ansatz as a linear combination of the usual 7-matrices in the Dirac representation 

/= (J _",) , f= , i- 1,2,3, (2.13) 

where a , a 2 ,a s are the Pauli matrices. In order to satisfy (2.2), we must transform these 
Dirac matrices of the vacuum into polar coordinates and multiply them by the factors T 
and yA, 

G t = T7 (2.14) 
G r = y/A (7 1 cos0 + 7 2 siwd cos 99 + 7 s sint? sin pj (2-15) 

= — (— 7 1 sin?? + 7 2 cos 1? cosy? + 7 3 cos?? sin 99^ (2-16) 

G^ = — — f-7 2 sin w + 7 3 cos . (2.17) 

r sin u V / 

This choice is convenient, because it greatly simplifies equations fl2.8D - (|2.1CD . Namely, the 
matrix p becomes independent of x and coincides with the usual "pseudo scalar" matrix 
7 5 in the Dirac representation, 



p = 7 5 = i 7 7 1 7 2 7 3 



1 

1 



As a consequence, the first and last summands in ( |2.9| ) vanish and thus 

B = ~ Tr(G m VjG n ) & G m G n 

~ Tr(G m V,G n ) (6l G n - b{ G m + CP g mn + ieK mnv 7 5 (?) . (2.18) 

Using Ricci's Lemma 

= 4V j5 ™ = V i Tr(G m G n ) = Tr((VjG m ) G n ) + Tr(G m (VjG n )) 

we conclude that the contributions of the first and second summands in the right bracket 
in ( 2.18j ) coincide and that the contribution of the third summand vanishes. Using the 
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anti-symmetry of the e-tensor, we can, in the contribution of the fourth summand, replace 
the covariant derivative by a partial derivative. This gives 

B = i Tr(G n Vj&) G n + -L e?™P Tr(G m £*,■<?„) 7 5 G P . (2.19) 

The second summand in fl2.19| ) is zero. Namely, the trace always vanishes if the tensor 
indices are all different, 

Tr(G m djG n ) = for m, j, n = t, r, i?, or (p and m^j^n^m ; (2.20) 

this can be verified directly using our special ansatz ( |2.14| )-( p.l7| ) for the Dirac matrices. 
In the first summand in fl2.19|) , we can use that VjG- 7 is a linear combination of the Dirac 
matrices G 3 , and thus 

Tr(G n VjG j ) G n = 4VjG j 

We conclude that 

B = ^jG j . (2.21) 

This form of B{ clS cl divergence of the Dirac matrices allows us to easily check that 
the Dirac operator is Hermitian with respect to the scalar product <. | .>; indeed 



In order to calculate the divergence ( |2.21|) , we first compute 



' dtU^G 1 ) = 



1 ff 
1 

Iff 



r T 



Iff 

= 0^ ( — 7 1 sin 2 $ + 7 2 sin # cos # cos 99 + 7 s sin ■& cos 1? sin </> 

r sin 77 V 

= f— 27 1 sini? cos i? + 7 2 (cos 2 $ — sin 2 i?) cos 02 + 7 s (cos 2 1? — sin 2 1?) sin 03 

r sin 77 v 

Ft Iff I G v ) = — ; — - (~7 2 cos ip - 7 s sin (p ) 

\n\ v r sin v \ / 



\g\ v rsini? 

and thus obtain 

i /2 T 



B = - [ ) G r — - (7 1 cos + 7 2 sini? cosy? + 7 s sin# sin</?) 

2 \r T ) r \ > 

= -n-A-hcr - - — G r 

r v ; 2 T 

We conclude that the Dirac operator has the form 

- ia i + a ' (4 + -r » - - 5 t) + iG 4 + iG 'h ■ (2 - 22) 
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3 The Dirac Equations 

In this section, we shall separate out angular momentum from the Dirac equation (J2T 
and reduce the problem to one on real 2-spinors. 

We first introduce some formulas involving Pauli matrices. These will be used in this 
section for the separation of the angular dependence, and then, in the next section, for 
the computation of the energy-momentum tensor needed in Einstein's equations. We 
introduce the following notation: 

<7 r ($,</?) = a 1 cos# + a 2 sint? cos <p + a 3 simd simp 
a (i?,<y?) = —a 1 simd + a 2 cosi9 cos ip + a 3 cost? sin ip 

(T (p ('&,(p) = (—(J 2 simp + a 3 cos 09) 

sin 17 

These matrices are orthogonal, 

Tr(cr r a®) = Tr(<7 r <7^) = Tr(</ o*) = 
and their square is a multiple of the identity, 



K) 2 

Furthermore, 

Tr (a* (8 v a r )) 
TrK (d#a r )) 



fi\2 



sin 2 # 



= 1 (3.1) 

2 sini? cost? (— cos ip simp + sin (p cos 99) = (3-2) 
cos $ 

(—sin 93 cos ip + cos 99 siny?) = . (3.3) 



sint? 



In analogy to the ansatz for the Dirac spinors in the hydrogen atom for zero angular 
momentum (see e.g. JTfl]), we write the wave functions in the form 



-iujt 



u\ e a 
a r u 2 e a 



1,2 



(3.4) 



where U\(r) and U2(r) are complex- valued functions, and the (e a ) a= i j 2 denote the standard 
basis e% = (1,0), e 2 = (0,1) of the two-component Pauli spinors. This ansatz is quite 
useful, because the Dirac equations for ^1 and ^2 are independent of each other, 



(T T 

-a r 



:VId r + - (VI- 1) 

r 

o 2i I a r 
+ wT7 +700 



2 T 



(3.5) 



where we have used (3.1). This allows us to view the Dirac equation as a two-component 
equation in m, u%. In order to simplify the radial dependence, we choose new functions 
$1 and $2 defined by 



$1 = r T 2 u\ , 



-ir T 2 u 2 



(3.6) 
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and rewrite the Dirac equation as 



1 
-1 



1 T 



1 

1 



Ad r + 



-1 

1 



m 



(3.7) 



This equation is real; thus we may assume that $ itself is real. On these two-component 
spinors, the scalar product <.|.> takes the form 



<$ I 



^ A 2 dr 



The "Dirac operator" in ( |3.7D is Hermitian with respect to this scalar product, and the 
normalization condition fl2.7|) for the wave functions transforms into 



I* 



2 T dr 



1 

47T 



Finally, we write the Dirac equation as the ODE 



A 



-1 

1 



1 

+ - 

r 



1 
-1 



m 



1 

1 



(3.8) 



(3.9) 



We remark that, instead of taking the ansatz (3.4) for the wave functions, we could 
just as well have put the matrix a r into the upper component, i.e. 



-iuit 



a' ui e a 
U2 e a 



1,2 



(3.10) 



This ansatz can be reduced to fl3,4| ) by changing the sign of the mass in the Dirac equation. 
To see this, we transform the wave function according to 

= 7 5 ^a 

Then since ^ is of the form ( |3.4| ) and satisfies the Dirac equation 
= 7 5 (G-m) 7 5 ^ a = -(G + m)^ a 



we can again simplify to the two-component equation (|3.7j) . We conclude that it also 
makes physical sense to look at the solutions of fl3.9| ) with negative m (and positive oS); 
they can be interpreted as solutions corresponding to the ansatz ( 3.10[) . 

In Appendix [A], we study the spinor dependence of general static, spherically symmetric 
solutions of the Einstein-Dirac equations, and we give a more systematic justification for 
the two ansatz' (|3.4|) and (3.10). 



4 Calculation of the Energy-Momentum Tensor 

We derive the form of the energy-momentum tensor by considering an arbitrary variation 
Sgij of the metric tensor. The variation of its inverse is given by Sg^ = —g lk g^ Sgki- In 
order to satisfy (^), we vary the Dirac matrices according to 

5GP = -\g^{5g kl )G l , 8G j = \ (5g jk ) G k . (4.1) 
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The energy-momentum tensor is obtained as the variation of the classical Dirac action 
(see e.g. @), _ 

S = J^(G-m)^ y/\g\d*x 

This action is real; therefore it suffices to consider the real part of the integrand. Since 
the wave function ^ solves the Dirac equation for the unperturbed Dirac operator G, we 
must only consider the variation of G; thus 

5S = fRe^(i{5G j ) 1 ^-+6B)^J\^\d 4 x . (4.2) 



We calculate the variation of the matrix B using (|4.l|) and (eTs|) , fl2.6|) , 

Re^5B^ = —Im5 ( Tr(G m V 7 -G n ) * G j G m G n *) 
16 V / 

= ^ 5 (Tr(G m VjG") Im (tt G j G m G n ^)) 
= _L ^ ( e r—p Tr(G m d,G n ) * 7 5 G P * 



16 

( = 0) ^e^^G™^) * 7 5 G P * 
1 



- e jmnp (*<7mfc) TV(G fc a,-G n ) * 7 5 G P * . (4.3) 



Notice that the trace in the last equation does not necessarily vanish. But we can use the 
fact that we are dealing with a spin singlet; this implies that the expectation value of the 
pseudovector 7 5 G p is zero, 

2 
0=1 

(this can be checked by a short explicit calculation). Thus we only have to consider the 



variation 5G J of the Dirac matrices in (|4.2| ). We substitute (4.1) into (4^2) and obtain for 
the spin singlet 



Thus the energy-momentum tensor has the form 

2 



a=l 



(4.4) 



We compute 

2 



Re ^* a iG*a t * a = 2wT|^| 2 = 2w T 2 r~ 2 |<£>| 2 

a=l 

2 _ _ 

Re ^¥^'G*<9 r * a = Re 2iT (r^ 1 Vf <9 r (r" 1 \/T $ 



a=l 
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Re £ iGr 9 t*o 

0=1 
2 

Re ^W^iCTd^o 



Rew^ ^ a G r ^ a = Re 2iuj 



Vat / o i 



o=l 
2 



-1 



$ = 



o=l 



(5.9) 



Re 2i^i (V 1 v^T $) ( _° . * j r (V 1 v^T $ 



Re 2^Irr 2 $ ( J q 1 ) 



-2w T 2 r~ 2 |$| 2 + 4T r~ 3 $i$ 2 + 2m T r~ 2 ($ 2 - $|) 

2 



Re iG * 9 ^*a = Re r" 1 £ * c 



a=l 



(3.1) 



Re J2 *a iGV W 



a=l 



(3.1) 



a=l 

-2r~ 3 T$i$2 

2 

Re r" 1 J]*; 

a=l 

-2r~ 3 T$i<J> 2 



il 
-il 



il 
-il 



All other combinations of the indices vanish because of ( |3.2[ ), ( |3 . 3[ ) , and the orthogonality 
of a r , a^, a^. We conclude that 



Tj = r~ 2 diag 2a; T 2 |<£>| 2 , -2a; T 2 \<S>\ 2 + 4Tr~ L $i$ 2 + 2mT ($f - 



-2T r" 1 $i$ 2 , -2T r" 1 $i$ 2 
As a check, we calculate the trace of Tj directly from ( |4.4j) , 

2 2 

r / = E Re d^Va) = ^ Re (^(iG^ +S)^) 



(4.5) 



a=l 



a=l 



0=1 



and we see that it agrees with taking the trace of (4.5) 



5 The Field Equations 

The Einstein Tensor Gj = Rj — 5j has the form (see e.g. [|l6| 



Cx 



1 


+ A + 


A' 






r 


1 


A 

+ ^2 " 


2AT' 




rT 


3 


A' 


AT 1 




27 ~ 


rT 



1/2 A rpff 



+ 
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with all other components vanishing. Thus, using (|4.5|), Einstein's field equations G 



-87T become 



16ttwT 2 |$| 2 = rA'-(l-A) 



(5.1) 



V 



-167rwT 2 |$| 2 + 327rr _1 T<£>i$2 + 16nmT (<S>j - <S>1) = 2rA — + (1 - A) (5.2) 



A 



T n A , T , 

r h r 

T 2 AT 



V 



T 



A' 



V 



2r I — — r — - + r — 



2A 



(5.3) 



as 



Using the notation a = $i, /3 = $ 2 > t ne equations (|3.9|), flO]), and (fy^j) can be written 



(5.4) 

(5.5) 
(5.6) 
(5.7) 



2rA 



A a 

Ap' 

r A' 
T' 



- a — (ujT + m) (3 
r 

{ujT — m) a — — (3 
r 

1 — A — WttujT 2 (a 2 + (3 2 ) 

A - 1 - 16itujT 2 (a 2 + /? 2 ) + 32vr - T a/3 + 16vr mT (a 2 - /3 2 



By direct computation (or, e.g. by using Mathematica) , we see that equation ( |5.3| ) is 
implied by this set of equations. It is sometimes useful to rewrite the Dirac equations 
(|5.4D and (5.5) in matrix notation as 



A& 



1/r —uoT — m 
ujT — m —1/r 



The normalization condition (2.7) takes the form 



$| 2 -U dr 
A 



1 

47T 



(5.8) 



If we assume regularity of the solution at r = 0, we obtain the following Taylor series 
expansions 

air + 0(r 3 ) 



a(r) 
0(r) 

A(r) 
T(r) 



1 



(ujTq - m) ai r 2 + 0(r 3 



1 wT 2 a 2 r 2 + C(r 3 ) 

3 



To 



m 
~6 



(AujT - 3m) T 2 a{ r 2 + 0(r 



(5.9) 
(5.10) 

(5.11) 
(5.12) 



where a%, To, u> and m are (real) parameters. We also require that our solutions have 
finite (ADM) mass; namely 



p := lim - (1 — Air)) < oo 



This implies that 



(5.13) 
(5.14) 



lim Air) = 1 

r^oo 

Finally, in order that the metric be asymptotically Minkowskian, we require that 

lim T(r) = 1 . (5.15) 
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6 Scaling of the Equations 



For the numerical computations, it is difficult to take into account condition ( 5.15] ) and 
the integral condition (^^). Therefore we find it convenient to temporarily replace these 
constraints by the weaker conditions 

lim T(r) < oo and / |<1>| 2 — == dr < oo (6.1) 
™ K ' JO ' ' y/A 

and instead set 

T = 1 and ra = ±1 . (6.2) 

This is justified by the following scaling argument which shows that there is a one-to-one 
correspondence between solutions satisfying (|Q|) , ( |BT2"D and solutions satisfying ( 5.15) ) and 

(PD- 

Thus, suppose we have a solution (a, /?, T, A) of ( |5.4| )-(|5?7D, ( p. 13 ) satisfying ( |6. 1[ ) , ( 
Then, with the parameters A and r defined by 



x T \ 5 



A = [4* J (o^)--* 



t = lim T(r) , 

r— >oo 

we introduce the new functions 

5(r) = y^a(Ar) 



/3(r) = ^/3(Ar) 

l(r) = A(Ar) 
f(r) = T~ l T(\r) 

A direct calculation shows that these functions satisfy the equations ( ^.4[ )-(pV^) with m, u> 
replaced by 

rh = Am , w = Xlot 



Moreover, it is easy to check that 

/ {dt 2 + p 2 )^=dr 
Jo \ A 



A 4vr 
lim T(r) = 1 

r— >oo 

and limr^oo | (1 — A(r)) < oo. We note that, as long as A,r are positive and non-zero, 
the transformation from the un-tilde variables to the tilde variables is one-to-one. 

Finally, we point out that for the physical interpretation, the conditions ( [5.8D and 
(|5.15|) must hold and hence only the scaled solutions can have physical significance. 



7 Existence and Properties of the Solutions 

Using the local Taylor expansion about r = 0, ( |5.9| )-( |5TT^ ), and setting T(0) = 1 and 
m = ±1, we construct initial data at r = 10~ 5 and then use the standard Mathematica 
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ODE solver. We shoot for numerical solutions by fixing a\ and varying oj. Using accuracy 
10 -8 , we found that solutions varied continuously with a±, (and to), indicating that this 
procedure is well-behaved. The solutions we find satisfy 



lim r 2 (a(r) 2 + P(r) 2 ) = , lim Air) = 1 

■ ->oo 



r 

lim — (1 — Ayr)) < oo , lim T(r) = r > , 

r^oo 2 r— >oo 

and, in addition, T and A stay positive for all r. In order to fulfill the important other two 
conditions ( |5.8| ),( ^15 ), we merely scale our variables in the manner described in section 



For any given a\ > 0, we found a unique ground state corresponding to the parameter 
value ujq, together with a countable number of distinct excited states corresponding to 
parameter values u n , n = 1,2,..., where Uo < u)\ < . . . < cj max (ai). For a\ = 0.02 
and m = 1, the scaled solutions for the ground state and for the first and second excited 
states are depicted in Figures |], ||| and |!|]. These solutions have the following general 
characteristics: The graphs of the functions T(r) have the same qualitative form for all 
values u n . The functions A(r), however, have their graphs changing with u) n ; indeed, for 
the nth excited state, the graph of A has precisely n + 1 minima. The tangent vector to 
the curve (a(r),/3(r)) for large r lies alternately in either the first or third quadrants. 

For larger values of ct\ (and m = 1), the solutions have a similar form, but the a — (3 
graphs of the ground state have self-intersections and are thus of a different shape. This 
is illustrated for three different values of a\ in Figure ^. We found that for the ground 
state, 

lim To(«i) = oo and lim To(«i) = 1 

Moreover, in every case, T is a monotone decreasing function of r. 

We also obtained solutions for m = — 1 by using similar methods. The results are 
qualitatively the same with the exception that, in this case, the tangent vector to the 
(a, (3) curve for large r lies alternatively in quadrants two and four, see Figures ||| and ^. 

The mass and the energy of the solutions we found have some interesting and surprising 
features, which we shall now describe. We consider the fixed n th excited state and, by 
varying a\ and adjusting u (for fixed m = ±1), we obtain a one-parameter family of 
solutions (parametrized e.g. by a\). After scaling, we find that solutions can only exist for 
a bounded range of m. For every value of m in this range, we obtain an at most countable 
number of solutions, which can be characterized by two physical parameters: the energy oj 



of the fermions, and the (ADM) mass ( 5.13 ). We find that we always have u < \m\, which 
means that the fermions are in a bound state. If we plot the binding energy u — \m\ vs. the 
mass m (i.e., the mass-energy spectrum), we find that this curve is a spiral which tends to 
a limiting configuration. This is shown in Figure ^ for the ground state (G), the first and 
second excited states (F, S), and the lowest and next excited states for negative mass (L, 
N). The spirals seem to have a self-similarity; this is illustrated in Figure |8] (we stopped 
the computation when the limitations of our numerics were reached). The (ADM) mass p 
can be viewed as the total energy of both the gravitational field and the fermions. Thus 
the quantity p — 2oj gives the energy of the gravitational field. As is shown in Figure |9[ 
it is always positive and, if plotted versus m, gives curves which again look like spirals. 
Furthermore, one can ask how much energy is gained (or needed) to form the singlet state. 

x We point out that we also plot the scaled variables in all the other figures. 



12 



For this, we must compare the (ADM) mass with the total rest mass 2\m\; i.e., the energy 
of two fermions at rest which are so much delocalized in space that their gravitational 
interaction becomes zero. This is plotted in Figure 10. For small m, the function p — 2\m\ 
is negative, which shows that one gains energy by forming the singlet state. (This gives 
a first hint that these states might be stable, because energy is needed to break up the 
binding.) For large values of m, however, (more precisely, shortly after the "turning point" 
of m; see Figure |H]), p — 2\m\ becomes positive. This indicates that the solutions should 
no longer be stable. This is indeed true and will be shown in the next section. 

We note that our scaling technique is essential for obtaining the mass spectrum - the 
unsealed variables do not have "spirals". 



8 Stability of the Solutions 

In this section, we shall consider the stability of our solutions under spherically symmet- 
ric perturbations. To this end, we consider the spherically symmetric, time-dependent 
Lorentzian metric of the form 



ds 2 



T- 2 (r,t)dt 2 - A-\r,t)dr 2 - r 2 {dfi 2 + sin 2 •& dip 2 



The time-dependent Dirac operator can again be calculated using ( 2.21D - One gets the 
following generalization of ( |2.22| ) (the dot denotes t-derivatives): 



G = G l 





'at 



4 A 



+ G r 



or r 



i_ r 

2 ~~T 



+ iCTde + iG*d« 



As in Section H|, we separate out the angular momentum by setting (cf. ( |3.4[ )) 



zi(r, t) e a 
i a r z 2 (r, t) e a 



a = 1,2 



This gives the two-component, time-dependent Dirac equation 
iTd t 

4 A 

( n i \ 

f Adr + 



1 
-1 



i m A i ■ x 
-T — + - T 
4 A 2 



1 

1 







-1 \ 1 



1 



111 







Observe that, in contrast to (|3.7|), this equation is complex, and thus we cannot assume 
here that the spinors are real. Calculation of the energy-momentum tensor and the Ein- 
stein tensor gives the equations 



1 A A' 
--+-+— 

T 2 A 
rA 

A 2AT' 



2iT 2 _ 
-87T — 5— Re (zi d t zi + z 2 d t z 2 ) 



-8tt Re 



iT 2 



(zi d r zi + z 2 d r z 2 ) + 



T^A-l 



(zi d t z 2 - z 2 d t zi) 



rT 



2TA? _ _ 
8tt — ^ — Re (z\ d r z 2 - z 2 d r z\) 



A' AT 1 A'T' 2AT AT" 3T 2 A 2 TAT T 2 A 
+ ^ * Tn + —-r- + 



2r rT 



T 2 



4A 2 



2A 



2A 



1tt — T Re(ziz 2 ) 
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Suppose that a solution (a(r), j3(r),A(r),T(r)) of the equations ( p.4|) -(p.7|) and (|5.q ), 
(|5.14 ), ( |5.15 ) is given for some values of the parameters m,u. Before making an ansatz for 
a time-dependent perturbation of this solution, it is convenient to introduce new spinor 
variables a(r,t), (3(r,t) by setting 

Zl (r,t) = e- iujt a(r,t) , z 2 (r,t) = e' iut /3{r,t) . (8.1) 



This has the advantage that we get the ansatz ( |3.4|) simply by assuming that a and j3 are 
time independent. We consider spherically symmetric, time-dependent perturbations of 
the form 

a(r,t) = a(r) + e(fli(r,t) +ia 2 (r,t)) (8.2) 
p(r,t) = P(r) + e(h(r,t) + ib 2 (r,t)) (8.3) 
A(r,t) = A(r) + eAi(r,t) (8.4) 
T(r,t) = T(r) + eT^t) , (8.5) 

where we look at the real and imaginary parts separately (oj and bj are real functions). 
Substituting into the Einstein-Dirac equations and only considering the first order terms 
in e gives (using Mathematica) a system of linear differential equations. If these equations 
admitted solutions with an exponentially growing time dependence, the original solution 
would be unstable. Therefore we separate out the time dependence by writing 

aj (r,t) = e Kt aj (r) , bfat) = e Kt b 3 {r) , j = 1,2, (8.6) 

A 1 (r,t) = e Kt A 1 (r) , Ti(r,t) = e Kt T^r) . (8.7) 

This gives the following system of five ordinary differential equations 

VAa[ = —-(m + uT)b 1 + KTb 2 --r^(--(m + ujT)(3]-uT 1 f3 (8.8) 
r 2 A \r J 

VAa'o = — - (m + uT)b 2 - kT h + k — T (3 - k — (3 (8.9) 
1 r 4A 2 

V~A b[ = — (m — wT) a\ — kT a 2 — ^ ^— (m — uT) a — —^j + to T\ a (8.10) 

V~Ab 2 = — (m — ujT)a 2 h kT a± — k —t Ta + k — a (8-H) 

2rAT[ = —^r- — Ti + ATi + (oi/3 + b x a) 

A r 

+16vrT 2 a (2ma 1 - 2uTa x + kTcl 2 ) - l§nT 2 (3 (2m&i + 2uTbi - nTb 2 ) 
-16?r Ti (^ujT 2 (a 2 + (3 2 ) - * T a[3 - 2mT (a 2 - (3 2 )^J 

+16vr ^j- (wT 2 (a 2 + p 2 )-^Ta(3-mT (a 2 - /3 2 )) (8.12) 

together with the algebraic equation 



Vat 

A 1 = 16vr- + 2u)b 2 )a + («»i + 2u)a 2 ) (5) . (8.13) 

KT 

The consistency of these equations (i.e., that the equation G\ = SirT^ is identically 
satisfied) was again checked with Mathematica. We want to show that there are no 
solutions for n > 0; this implies stability. 
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The above equations come with initial conditions at r = and additional constraints, 
which we will now describe. A Taylor expansion about r = gives, similar to 
the following expansions near r = 

ai(r) = a w r + 0{r 2 ) , a 2 (r) = a 20 r + 0(r 2 ) (8.14) 

6i(r) = 0(r 2 ) , 6 2 (r) = 0(r 2 ) (8.15) 

A x {r) = 0{r 2 ) , Ti(r) = T 10 + 0(r 2 ) . (8.16) 

We have three parameters a±o, a 2 o and Tio to characterize the solutions. Since the metric 
must be asymptotically Minkowskian, we demand moreover that 

lim AUr) = (8.17) 

1 — »oo 

lim Ti(r) = . (8.18) 

r — >oo 

Furthermore, the wave functions must be normalized, which means that (cf. (|3,8])), 

2 T(r t) 1 

{a 2 (r,t) + ^t))^L )d r = - , (8.19) 

for all t. This time-dependent normalization condition appears to make the analysis very 
complicated. It turns out, however, that we do not have to consider it at all, because, for 
the perturbation (jsl|), (j87?|), it is automatically satisfied. 

To see this, note that, as a consequence of the current conservation (2.6), the nor- 
malization integral ( 8.19j ) is actually time-independent. But, in the limit t — > — oo, the 
time-dependent solution (|8.2[) -( ^5|) of the Einstein-Dirac equations goes over into the 
static, unperturbed solution (a, (3, A,T), and thus ( |8.19| ) holds in this limit. It follows 
that fl3.19|) holds for all t. In other words, the linear contribution in e to ( |8.19| ) vanishes 
as a consequence of the linearized Einstein-Dirac equations (|8.8|)-( |S7L2| ), ( |3.13| ). We must 
only make sure that the integral ( 8.19 ) is finite for all t. 

These conditions can be simplified. To see this, we first consider the infinitesimal time 
reparametrization 

k T(r = 0) 

This transformation does not destroy the ansatz ( |8.2| )-( |8?7|) ; it only changes the functions 
a 2 , 6 2 and T\ according to 

Tl(r) ^ Tl(r ) _ ^(0)^1 

a 2 (r) - a 2 (r) - ^ ^ a(r) (8.20) 

Hr) ^ Ur ) - ^^P(r) . (8.21) 

Thus we can arrange that T\ vanishes at the origin, 

Ti(r) = 0{r) 
provided that ( |8.18| ) is replaced by the weaker condition 

lim T\(r) = /j, for some fi, < /x < oo 
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This makes the numerics easier, because we now have only two free parameters axo, C120 at 
the origin to characterize the solution. Furthermore, using the linearity of the equations, 
we can multiply every solution by a (non-zero) arbitrary real number. This allows us to 
fix one of the parameters (e.g. by setting 020 = 1); an d thus we end up with only one free 
parameter. 

Our strategy is to show that, for any k > 0, there are no solutions for which the 
normalization integral ( |8.19 ) is finite; this will imply stability. In order to explain the 



technique and the difficulty for the numerics, we consider Figure 11, where typical plots 
of (01,61), (02,62) for a small value of k, and the ground state solution of Figure [l], are 
shown. According to ( |8.14|) , (|8.15 ), both the (a±,bx) and the (02,62) curves start at the 



origin. We want to show that at least one of these curves stays bounded away from the 
origin as r —* 00, no matter how we choose k and the initial values aio, 020- This will 
imply that the integral ( |3.19| ) is unbounded. Figure 11 is interesting because it almost 



looks as if this happened: the (02, 62) curve looks like the a — /3-plot of the ground state, 
whereas (ax, bi) is almost like the a — /3-plot of the first excited state. The (ax, 61) curve 
is not quite closed, however, and it is difficult to decide whether this is just an artifact 
of the numerics or whether it actually means that there is no normalizable solution. The 
numerics are especially delicate because, for small values of k, the values of (02,62) are 



much larger than (ai,b\) (notice the different scales on the plots of Figure 11). In order 
to improve the accuracy of the numerics, it is useful to note that, numerically, 02 and 62 
are, to a very good approximation, multiples of the unperturbed wave functions a, (3, i.e. 
02 ~ /ia, 62 ~ for some real constant [i. We can eliminate this dominant contribution 
to ci2, 62 by introducing new variables 

02 = a 2 - n a , 6 2 = 6 2 - fi /3 

and rewriting our ODEs in the functions (ax,bx,a 2 ,b2,Ax,Tx); this gives a system of five 



inhomogeneous, linear equations. From these, we obtain the plots in Figure 12. The 



(ax, 6i)-curve is similar to that of Figure 11; the (0,2, 62)-plot, however, gives a much more 
detailed view of the imaginary part of a(r,t), (3(r,t) (notice again the different scales on 
the plot). 

Next we describe our method to determine the initial data at r = for the solutions 
(i.e. the parameters aio, O20 5 for the original equations (|3.8|)-( |S7r2"| )). It turns out that the 
integral ( |8.19| ) is only finite for all t if both (ai,&i) and (02,62) tend to zero for large r; 
indeed, from the numerics we see that (ax, 61) and (02, 62) cannot have infinite oscillations. 
It is an efficient technique to fix the initial values by trying to satisfy the first part of these 
conditions 

lim(ai(r), bx(r)) = (0,0) . (8.22) 

r — >oo 

This can be done by varying the initial values in such a way that 

(ax(R) 2 + a 2 (Rf) -» min 

where R is the value of r where we stop the numerics (R must be chosen sufficiently large). 
Using the linearity of our ODE's, this minimizing condition leads to simple algebraic 
relations between the initial data and the values of ax(R), b±(R) for a fundamental set of 
solutions, from which the initial values can be computed numerically. According to Figures 
pd| and |l^, we already know qualitatively how the resulting (ax, &i)-plot is supposed to look 
like. This is helpful for checking the numerical results and for finding the best values for R. 
(More precisely, the best value for R is, roughly speaking, the value where the (ax, 6i)-plot 
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starts to look like a closed curve. If R is chosen much larger than this, the numerical 
inaccuracies pile up in a such way that our method of computing the initial data from 
ai(R) and bi(R) is no longer reliable.) 

This procedure can be carried out for different values of n, and gives the plots of 
Figures 13 to [0^. For very small values of k, the plots look like those in Figure |l3], and 
one sees that (pair), 62(f)) does not tend to zero for large r. The shape of the (S2, £>2)-plot 
does not change over many magnitudes in k (see Figures 12, |l4|, |D^), showing too that 
our numerics are well-behaved. For k ~ 0.02, the form of the plots changes drastically 
(see Figures [H], 17, |l^). One sees that (02,62) still do not go to the origin and that it 
becomes impossible to satisfy even condition (|8.22|) . If k is further increased, both (ai, b\) 
and (02,62) g° to infinity for large r, no matter how the initial data is chosen. From this, 
we conclude that the ground state solution of Figure [l] is linearly stable. 

Our method also applies to the excited states. For the first excited state and the 
lowest negative-mass state (i.e., for the solutions of the Figures Q and |5|), the solutions of 
the linearized equations for small k are plotted in Figures [l9| and |2^. The (02, ^-curves 
are again of similar shape as the corresponding a — /3-plot; the (a±, 6i)-curve resembles the 
a — /3-plot for the next higher excited state (i.e., in Figure [ji] for the second excited state, 
and in Figure 20 for the first excited negative-mass state). It is again useful to introduce 
the variables 0,2, 62- A numerical analysis of the equations for different values of n shows 
that these solutions are also stable. 

It might seem a bit surprising that even the excited states are linearly stable. Actually, 
this can already be understood qualitatively from Figure 1£ and the following heuristic 
argument. Thus, if the solution (for small coupling) were unstable, then the solution 
(a,j,bj,Ai,Ti) of the linearized equations (^q )-( 8T3|) would, to first order, describe the 



decay of the bound state. Therefore the (aj, 6j)-plots give us information into which state 
the wave functions tend to decay. The (a2,i>2)-plot is not interesting in this respect, be- 
cause it looks like the a — /3-plot and only yields information about a time-dependent phase 
transformation of the wave functions. The (01, 6i)-plot, however, is helpful. According 



to Figure 19, the wave function tends to decay into the second excited state. This is 
surprising; one might instead have expected the tendency to decay into the ground state. 
Since the energy of the second excited state is higher than that of the first excited state, it 
would seem physically reasonable that this decay cannot happen spontaneously; this gives 
a simple explanation for stability. 

We point out that these stability results are only valid for weak coupling, i.e., for small 
(scaled) mass m. If m comes into the region where the mass spectrum of Figure ^ starts to 
have the form of a spiral curve, the numerical behavior of the linearized equations becomes 
much more difficult to analyze and does no longer allow simple conclusions. It is thus very 
helpful to study the stability in this regime via a different method, which we will now 
describe. 

The existence of the spiral curve (c.f. Figure [0]), enables us to obtain information 



regarding the stability properties of these solutions, using Conley Index theory (see |14 



Part IV]). Indeed, as the figure shows, if m > mi, there are no solutions while at m = mi, 
the solution Pi suddenly appears. This solution is "degenerate", and has Conley index 
(the homotopy type of a one-point, pointed space). For m < mi, the solution P\ 
bifurcates into two solutions Qi and Qi- The solution Qi, being a "continuation" of 
the stable solution Qq (m is near 0), must also be stable; this follows from Conley's 
Continuation Theorem |14|, Thm. 23.31]. In fact, since Qq being stable implies that the 
Conley index of Qq is £ , (the pointed zero sphere), the continuation theorem implies that 
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the Conley index of Q\ must also be E , and thus Q\ is also a stable solution. (Moreover, 
the same argument shows that all points on the curve between and Pi correspond to 
stable solutions; this can also be checked numerically). Since the Conley indices of Q\ and 
Q2 must "cancel" (i.e., the index of Q\ and Q2 taken together must be 0), this implies 
that Q2 cannot be stable, and in fact, the index of Q2 is E 1 , the pointed 1-sphere. (In 
fact, all points on the curve between Pi and P2 correspond to unstable solutions; one can 
also check this numerically.) Similar reasoning can be applied to solutions near P2, P3, 
. . . , and so on. (We remark that only m, and not u>, can serve as a bifurcation parameter; 
we show this in Appendix |B|.) 

These general Conley Index techniques also enable us to show that for each n, these 
spiral curves must tend, as a\ — > 00, to a limiting configuration r n , which is either a 
single point, or is a "limit cycle" S l ; i.e. a topological circle (we assume, as the numerics 
indicate, that the curve "spirals inwards"). In fact, were this not the case, then for each 
value of m between and mi, the corresponding solution set would form an "isolated 
invariant set" (c.f. fl4|| ), and so, again by Conley's Continuation Theorem, their Conley 
indices would all be the same. However, for m near 0, the index of the isolated invariant 
set is E°, while for m near mi, the corresponding isolated invariant set has index 0. Since 
these two indices are different, we have a contradiction. 

It follows from this last result that for a point (m, m — Cc?) on r?7, there are an infinite 
number of solutions with m = m, as well as an infinite number of solutions with m — u = 
rh—uj. For parameter points not meeting T n , there are at most a finite number of solutions. 



A Appendix: Justification of the Ansatz for the Spinor De- 
pendence 

In this section, we consider the general form of the spinors in static, spherically symmetric 
systems and derive the Einstein-Dirac equations for these systems. This analysis will 
clarify the ansatz' (|3.4j) and ( |3.10D for the wave functions, which was made in Section ||| 
without a detailed explanation. 

The Dirac wave functions $1, ^2 of a general two-fermion system can be written in 
the form 

V a {x,t) = A{x,t)e a , a = 1,2 . (A.l) 

where A® = VP" is a (4 x 2)-matrix and where (e a ) again denotes the standard basis of the 
two-component Pauli spinors. The system being static means that the time dependence 
of A has the form of a plane wave, 

A(x,t) = e~ iujt A(x) 

The simplest way to characterize the spherical symmetry of the spinors is to demand that 
the angular dependence of A is described only by the submatrices 1 and a r ; i.e. 

ui(r) 1 + v 2 (r) a r \ . , 

v 3 (r) 1 + U4 (r) a r ) { > 

with complex functions Vi, . . . , V4. This form of A can be derived if one requires that the 
total angular momentum is zero and that all the expectation values <0> = J2a=i ^aO^a 
of the spin matrices are spherically symmetric. In a simplified argument, this form can be 
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understood directly from the fact that the presence of any matrices a , a v would destroy 
the radial symmetry in ([A. 2 ). 

The ansatz ( [0] ), ( JO) ) for the wave functions is a linear combination of Q3.4| ) and 
( ft.iop . From this, we immediately obtain the corresponding Dirac equations. Namely, the 
complex two-spinors <3? = ($i,$2) and S = (Si, S2) with 



$1 
$ 2 



rT 2 vi 
_ 1 

—ir T 2 V4 



rT 2 V2 
_ 1 

— ir T 2 y 3 



satisfy, (in analogy to ( p^9|) ), the equations 
A* ~' ° 

A S' 




1 



-1 





1 



- 1 

r 



-1 




(A.3) 
(A.4) 



In Section ||, we argued that the reality of the coefficients in ( |3.9| ) allows us to choose real 
spinors. This procedure simplified the Dirac equations considerably, but it might be too 
restrictive to describe the general solution of the Einstein-Dirac equations. In order to 
analyze the situation more carefully, we first note that the function 



F(r) 




$(r) 



is independent of r, as is obvious from ( A.3 ). The boundary conditions at the origin, 
3>i(0) = = ^(O), imply that F must vanish identically, and thus the product ^1^2 is 
real. This means that $1 and $2 must be real up to a common phase factor, i.e. 



$i(r) 



fl(r) 



$ 2 (r) 



/ 2 (r) 



(A.5) 



with real functions f\ , / 2 . Again as a consequence of the Dirac equation QA.3| ) , the phase 
a is independent of r. A similar argument applies to S and gives 



=1 



9i(r) 



-<2 



92{r) 



(A.6) 



with real functions g\ , 52 and (3 G IR. 

We shall now compute the energy-momentum tensor. First of all, the spherical sym- 
metry of the spinors implies that the off-diagonal components T$, T* T%, T r n , Tz vanish 



and that T$ = T£. Thus we must only consider T* and the diagonal components Tf , T£, 
Tt of the energy-momentum tensor. As a first step, we verify that the contribution (HT 



of the variation of B vanishes: According to ( [2.20 ), the trace in ( |4.3| ) is zero if 5g m k is 
diagonal. Thus we must only consider if, and we may therefore assume that the indices 
m, k are either m = t, k = r or m = r, k = t. Furthermore, the spherical symmetry ( |A.2| ) 
implies that the expectation of the pseudovector 7 s G p only has a component in the time 
and radial directions, 



7 5 G P fy a = for p = <p or p = •& 



a=l 
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We conclude that both indices m and p in (|4.3| ) must be equal to either r or t, and the 
antisymmetry of the e-tensor implies that the remaining indices j, n must coincide either 
with § or tp. Thus we must only consider the trace in (4.3) for the combination 



Tr (G k (d#G v -d v G<>] 



(A.7) 



But, according to ( 2.16 ), ( 2.17| ) and ( 2.11 ), we have d^G v = d v G^, so that ( A.7 ) vanishes. 
We conclude that the energy- momentum tensor is again given by (4.4). In order to 



compute the trace in (4.4), it is useful to first notice that, if we write the Dirac matrices 
as (2 x 2) block matrices, then G l is diagonal with entries proportional to the identity, 
whereas G r is off-diagonal with submatrices which are multiples of a T . This implies that 
the mixed contribution (i.e., the contribution proportional to <&E or H<£) to Tf and TJT 
vanish. Using the explicit form of G^ together with (3.1), we conclude that the mixed 
contribution also vanishes in T$. Thus the energy-momentum tensor of the system is 
simply the sum of the energy-momentum tensors corresponding to the spinors $ and S. 
As a consequence, the constant phase factors in ( A.5D , ( |A.6|) are irrelevant; we can without 
loss of generality assume that $ and E are real. Using (|4.5| ), we end up with the formulas 



rp! 

3 



T][$] + Tj\ 



with 



(A.8) 
(A.9) 



7][$] 



diag (2wT 2 |$| 2 , —2u T 2 \&\ 2 + AT r" 1 + 2mT($ 2 -$ 2 ,), 

-2Tr~ l $!$ 2 , -2Tr~ l $ a S a ) (A.10) 
diag fou T 2 |H| 2 , -2lo T 2 \E\ 2 - AT r" 1 E{E 2 + 2m T (H 2 - S|), 

(A.ll) 



2Tr- 1 E 1 E 2 , 2Tr- 1 E 1 E 2 



Thus the Einstein-Dirac equations take the form flA.3|) ; (|A.4|) and 



-8tt (T][$]+T|[S] 



with Tj given by ( A. 10 ), ( A.ll ). This is a generalization of the system (|5.4j)-(5.7). It is 
quite remarkable that the energy-momentum tensor is just the sum of Tj[$] and 
Similar to (^]B), the normalization condition for the wave functions takes the form 



(1*1 



2 _|_ p|2 



T 



A 



dr 



1 

Att 



(A.12) 



We now qualitatively describe how the solutions of this generalized system can be 
constructed and how we recover the solutions of Section [7]. The scaling technique of 
Section ||| can again be applied and consequently we can omit the conditions T(oo) = 1 
and ( |A.12| ) if we instead set T(0) = 1 = m. Then the solutions are characterized by 
the three parameters u, $i(0), 2^(0) (notice that a Taylor expansion around r = yields, 
in analogy to ( |5.10| ), the constraints ^(O) = = E^(0)). Compared to the situation in 
Section ^, we thus have one additional continuous parameter to describe the solution. At 
first sight, this might seem to imply that we can now construct, for given u>, a continuous 
one-parameter family of solutions. Then our ansatz' ( |3.4|) and ( |3.10j ) would just correspond 
to two special points of this continuum of solutions, and it would become unsatisfying that 
we just picked these two special solutions for the discussion of the mass spectrum and the 
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stability. However, the additional free parameter is illusory due to the fact that we also 
have one additional constraint at infinity. Namely, we must (for given lo) adjust ^(O) and 
H' 2 (0) in such a way that both the ($1, $2) and the (Hi, H2) curves tend to zero for large 
r. For generic u>, these conditions will only be satisfied for a discrete set of initial values 
($i(0),H' 2 (0)). The choices ($i(0) = 0, H' 2 (0)) and (*i(0), H' 2 (0) = 0) are both allowable; 
they correspond to the solutions constructed in Section [7|. After scaling, this shows that 
for generic m, the Einstein-Dirac equations only admit a discrete number of solutions. 

We note that these considerations do not rule out the possibility that the general 
ansatz for the spinors might lead to some additional solutions. We did not study the 
general equations systematically and can only say qualitatively that it seems difficult 
to arrange that there are simultaneous normalizable solutions of (A. 3) and (A. 4). The 
existence of solutions of this type, however, remains an open question. 



B Appendix: Justification of Using m as the Bifurcation 
Parameter 

We show in this section, first that lo is unsuitable as a bifurcation parameter, and second 
that m can serve as a bifurcation parameter. In Conley Index theory, a parameter can only 
be used as a bifurcation parameter if it remains well-defined and fixed when perturbations 
of the solutions are considered. The basic reason why m and lo must be treated differently 
can already be understood from the general form of the Dirac equation in (|l.l|). The mass 
m enters as an a-priori given parameter into the Dirac equation, whereas the energy lo 
of the fermion is only determined by the solution vl/. This means that if we consider a 
variation of a solution, m can be considered as a fixed parameter, while lo will in general 
change. If the perturbation of the solution is not static, the energy of the fermion will in 
general become time-dependent, and lo will no longer be a well-defined parameter. 

Although being correct in principle, this argument is too simple and not fully con- 
vincing. First of all, the situation becomes more complicated by our scaling technique, 
which also changes m and thus makes it impossible to consider the mass as a fixed pa- 
rameter throughout. Furthermore, lo is uniquely determined by the solutions (a, f3, A, T) 
of ( |5.4| )-(pV7|). It enters as a parameter into the linearized equations in a similar way as m 
does, and it is not obvious from these equations why the two parameters m and lo should 
play such different roles for stability considerations. Therefore we will show in detail that 
solutions of the linearized equations do not determine lo, whereas m is still a well-defined 
parameter. 

In order to show that lo is not well-defined, we generalize the ansatz of Secti on || in the 
way that we also allow lo to be time-dependent. In analogy to fl3.2j)-( |8l| ) and (|8.6[) ,(8.7), 
we consider a perturbation of lo of the form 

to(t) = LO + €LO ie Kt . (B.l) 

Since lo represents a frequency, i.e. the time-derivative of a phase, the correct generalization 
of equation flS.l|) is to replace the phase factor e~ lujt by 

exp (— i J io(s) ds 

The ansatz for aj, bj, A\, T\ then remains the same as before, given by (^.2|)-( |8^5| ) and 
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08.6D , (|8.7|) . Thus the spinors z±, z 2 are given by 



zx(r, t) 
z 2 (r, t) 



■*fi u W da [a(r) + s{a 1 (r,t) + ia 2 (r,t)} 
^ w(,,d, ^(r) + e(6i (r, t)+i 62 (r,t)] 



Into these equations, we substitute ( |B.l[) and consider only the first-order terms in e. This 



gives 



Zl (r,t) = e~ iwt (l - ie J 001 e KS ds^j[a(r) + e(ai(r,t) + ia 2 {r,t)} 



-iujt 



-iujt 



air) 



ie — (e Kt - 1) a(r) + e(ai(r,t) + i a 2 (r,t) 

K 



(1 + ie — ) a(r) + e(ai(r, t) + i 02(7, t) — ie — e Kt a(r) 



-i(wt-e ^) 



a(r) + e(ai(r,t) + i a 2 (r,t) — ie — e Kt a(r) 

K 



with a similar expression for z 2 (r,t). This looks quite similar to the original ansatz (8.1) 
and ( |8,2|) -( |8T5|) except for two differences; namely, there is here an additional constant 
phase factor exp(ie — ), together with the term 



ie — e Kt a(r) 

K 



(B.2) 



The phase factor exp(ie ^) plays no role in our analysis, since it falls out of all the 
equations (notice that the energy-momentum tensor (4.4) contains only terms of the form 
What is interesting about the term (|B.2| ), however, is that its time-dependence 
is again of the form e Kt . It is thus consistent with our ansatz ( |8.6| ), and corresponds to 
the transformation 



and similarly 



a 2 {r) -> a 2 (r) 



b 2 {r) -» b 2 (r) 



e — air) 

K 



(B.3) 



(B.4) 



Thus the more general ansatz for a,j, bj,Ai,T\ whereby co is replaced by (p3. 1|) is the same 



as the original ansatz for these quantities if we transform a 2 and b 2 according to (B^) and 
(B.4). Conversely, we may obtain an arbitrary time dependence in u, of the form ( |B.1| ), 
merely by transforming a 2 and b 2 according to (|B.3| ) and (|B.4| ) . This means that a solution 
of the linearized equations only determines u up to linear time-dependent perturbations 
of the form (BT). Thus uj is not a well-defined parameter^]. 

For the parameter m, however, the situation is completely differently. Namely if ^ is 
a solution of the time-dependent Dirac equation = m^, then we see from fl2.7| ) that 
(Gty I G^>) = m 2 . But as we noted earlier, (after (|8.1S| )), this relation is also valid for the 
linearized equations. That is, m is a well-defined parameter for the linearized equations. 



2 It is interes ting t o not ice that the contribution to ((12,62) proportional to (a, (3) which occurs in the 
transformation (B.3), (B.4) played an important role in our numerics. Namely, we saw in Section H that 



the (a2,62)-plot looks very similar to the (a, /3)-plot (see Figure [Tl]), which shows that this contribution 
is actually dominant for small k. It caused problem s in t he n umerics and forced us to introduce the 
new variables 0,2, 62 (see Figure [12[) . According to (B.3),(p3^4[), we can now understand the do mina nt 



contribution to (02, 62) in Figure [fT|as describing a time-dependent perturbation of uo of the form (B.l ) 
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Figure 1: The ground state for m = 1, a.\ = 0.02; scaled parameter values: m = 0.5340, 
u = 0.4994 




Figure 2: The first excited state for m = 1, ol\ = 0.02; scaled parameter values: m = 
0.7779, w = 0.7326 



beta 




Figure 3: The second excited state for m = 1, a± = 0.02; scaled parameter values: 
m = 0.9616, uj = 0.9080 
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Figure 4: a-/3-plot for the ground states at parameter values m = 1 and a.\ 
0.31,0.35,0.45 




Figure 5: The lowest state for m = —1, ol\ = 0.02; scaled parameter values: m = —0.7567, 
to = 0.6302 



beta 




Figure 6: The next excited state for m = —1, ol\ = 0.02; scaled parameter values: m = 
-0.9742, uj = 0.8391 
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Figure 8: Binding Energy \m\ — co of the Ground State, Detailed Pictures 
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Figure 15: Perturbation of the ground state for k = 0.01 
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Figure 16: Perturbation of the ground state for k = 0.02 
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Figure 17: Perturbation of the ground state for n = 0.023 
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Figure 18: Perturbation of the ground state for k = 0.03 





Figure 19: Perturbation of the first excited state for k = 10 6 
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